def warn(*args, **kwargs):
pass
import warnings
warnings.warn = warn
import os
import numpy as np
import pandas as pd
import scanpy as sc
import scvi
import matplotlib.pyplot as plt
import seaborn as sns
sc.settings.verbosity = 3
sc.logging.print_header()
sc.settings.set_figure_params(dpi=100, facecolor='white')
Global seed set to 0
scanpy==1.9.1 anndata==0.8.0 umap==0.5.3 numpy==1.22.3 scipy==1.9.3 pandas==1.5.2 scikit-learn==1.1.3 statsmodels==0.13.5 python-igraph==0.10.2 pynndescent==0.5.8
adata = sc.read_h5ad('./data/retinal.h5ad')
CellAssignによるCell typeの自動推定。
Zhang, A.W., O’Flanagan, C., Chavez, E.A. et al. Probabilistic cell-type assignment of single-cell RNA-seq for tumor microenvironment profiling. Nat Methods 16, 1007–1015 (2019).
CellAssignは、事前知識としてCell typeの種類と、そのCell typeを識別するマーカー遺伝子の組み合わせの情報を使って、それぞれの細胞がいずれのCell typeに割り当てられるのか(あるいは"unassigned"か)の確率分布を計算してくれるツール。
観測された値が、cell type, library size, batchの影響を受けると仮定。
入力は遺伝子のカウントと、マーカー遺伝子・cell typeの対応のテーブル。発現量の閾値などの設定は必要ない。マーカーとcell typeの対応が0/1で表現されていればオーケー。
基本的には、scVIと似たような負の二項分布のモデリング。バッチIDやライブラリサイズを取り込んで、細胞ごとにcell typeを表す潜在変数の確率分布をEMアルゴリズムで最尤推定する。オリジナルのR版はtensorflowで最適化しているが、scvi-toolsに取り込んで実装したバージョンではpytorchを使って最適化。
マーカー遺伝子の情報は、エキスパートの知識で与えるか、CellMarker、PanglaoDBなど、マーカー遺伝子のデータベースからとってきて用意する。
ここでは、CellMakerからダウンロードしたマウスのマーカー遺伝子全情報をまとめたエクセルファイルをダウンロードして、該当する組織(眼球・脳・胚)とアノテーションしたいCell typeに対応する遺伝子のリストを抽出し、0/1のテーブルにまとめている。エクセルファイルのサイズがでかいので講習では割愛。すでに処理済みのeye_markers.csvのテーブルを配布している。
%%script true
target_labels = ['Amacrine cell',
'Cone photoreceptor cell',
'Horizontal cell',
'Late activated neural stem cell',
'Late neuroblast',
'Neural progenitor cell',
'Neuroblast',
'Photoreceptor cell',
'Retinal ganglion cell',
'Retinal progenitor cell']
all_markers = pd.read_excel('./models/marker_genes/Cell_marker_Mouse.xlsx')
eye_markers = all_markers[(all_markers['tissue_class'] == 'Eye') | \
(all_markers['tissue_class'] == 'Brain') | \
(all_markers['tissue_class'] == 'Embryo')]
eye_markers = eye_markers[['Symbol', 'cell_name']]
eye_markers = eye_markers[eye_markers['cell_name'].isin(target_labels)]
eye_markers['dummy'] = [1]*len(eye_markers)
eye_markers = eye_markers.dropna()
eye_markers = eye_markers.drop_duplicates(keep='first')
eye_markers = eye_markers.pivot(index='Symbol', columns='cell_name', values='dummy')
eye_markers = eye_markers.fillna(0).astype(int)
eye_markers.to_csv('./models/marker_genes/eye_markers.csv')
eye_markers = pd.read_csv('./models/marker_genes/eye_markers.csv', index_col=0)
以下のように、マーカー遺伝子対Cell typeの対応関係が0/1で表現されている。
eye_markers.head()
| Amacrine cell | Bipolar cell | Cone cell | Cone photoreceptor cell | Horizontal cell | Late activated neural stem cell | Late neuroblast | Müller glial cell | Neural progenitor cell | Neural stem cell | Neuroblast | Neuron | Photoreceptor cell | Retinal ganglion cell | Retinal progenitor cell | Rod photoreceptor cell | Rod photoreceptor precursor cell | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Symbol | |||||||||||||||||
| 4932438A13Rik | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| 6330403K07Rik | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Aars2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Abca4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 |
| Abca5 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
CellAssignはライブラリサイズの情報が計算に必要。ここでは細胞に割り当てられたトータルリード数の、全細胞平均に対する比率で計算することにする。前処理段階で計算した 'total_counts' を流用しているが、カウントデータを使って足し算を計算しなおしてもいい。対数変換した値で足し算しちゃうとスケールが変わっておかしなことになっちゃうので注意。
adata.obs['size_factor'] = adata.obs['total_counts'].values / adata.obs['total_counts'].values.mean()
マーカー遺伝子のうち、高発現変動遺伝子とオーバーラップするものだけ抽出。(計算速度のため。時間に余裕があれば発現変動遺伝子に限定せず全部のマーカー遺伝子を使ってもいい)
該当のマーカー遺伝子の情報だけ取り出した新しいanndataオブジェクトを用意する。
eye_markers = eye_markers[np.isin(eye_markers.index, adata.var[adata.var['highly_variable'] == True].index)]
bdata = adata[:, eye_markers.index].copy()
モデルのセットアップ。ライブラリサイズやバッチに対応したラベルを指定して、カウントデータが格納されているレイヤーをセットする。
scvi.external.CellAssign.setup_anndata(bdata,
size_factor_key='size_factor',
batch_key='batch',
layer='counts')
WARNING:jax._src.lib.xla_bridge:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
学習の実行。やはり時間がかかるので講習では学習済みのモデルを配布。
%%script true
cellassign_model = scvi.external.CellAssign(bdata, eye_markers)
cellassign_model.train()
%%script true
cellassign_model.save('./models/cellassign/', overwrite=True)
cellassign_model = scvi.external.CellAssign.load('./models/cellassign', bdata)
INFO File ./models/cellassign/model.pt already downloaded
細胞ごとのcell type確率の予測。
predictions = cellassign_model.predict()
predictions
| Amacrine cell | Bipolar cell | Cone cell | Cone photoreceptor cell | Horizontal cell | Late activated neural stem cell | Late neuroblast | Müller glial cell | Neural progenitor cell | Neural stem cell | Neuroblast | Neuron | Photoreceptor cell | Retinal ganglion cell | Retinal progenitor cell | Rod photoreceptor cell | Rod photoreceptor precursor cell | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 3.376599e-04 | 5.895629e-06 | 9.812665e-01 | 1.480903e-03 | 3.604459e-05 | 5.409688e-04 | 9.383034e-04 | 5.445538e-04 | 1.180385e-04 | 2.714216e-15 | 1.335013e-24 | 6.517411e-106 | 1.963789e-03 | 5.391997e-03 | 6.316957e-03 | 1.042139e-03 | 1.623581e-05 |
| 1 | 8.178778e-04 | 3.013952e-03 | 9.781281e-01 | 2.100749e-03 | 5.616672e-05 | 1.265442e-04 | 1.375039e-03 | 9.361648e-04 | 4.931543e-05 | 4.166674e-16 | 1.365678e-21 | 1.209643e-83 | 3.419927e-03 | 2.560188e-03 | 6.296753e-03 | 1.038806e-03 | 8.046374e-05 |
| 2 | 2.660053e-16 | 4.737028e-16 | 2.269205e-11 | 6.727371e-19 | 4.591848e-26 | 1.872102e-33 | 5.113711e-22 | 1.136622e-21 | 5.638754e-21 | 3.098262e-28 | 2.559535e-51 | 1.390883e-139 | 6.139165e-15 | 1.000000e+00 | 1.460813e-13 | 2.409974e-14 | 3.162411e-22 |
| 3 | 4.385756e-07 | 2.792256e-09 | 5.797536e-04 | 1.120042e-06 | 1.672435e-07 | 5.477194e-08 | 5.914494e-11 | 1.707186e-07 | 2.253600e-07 | 1.501606e-09 | 6.281428e-21 | 6.914016e-82 | 1.692987e-06 | 9.994120e-01 | 3.732196e-06 | 6.157184e-07 | 3.138582e-08 |
| 4 | 1.568833e-05 | 8.424272e-05 | 9.923665e-01 | 7.462969e-07 | 4.404270e-17 | 1.324441e-12 | 6.180864e-07 | 2.804698e-12 | 4.296175e-07 | 1.798462e-28 | 4.651498e-44 | 1.642679e-157 | 8.939292e-05 | 9.253567e-13 | 6.388414e-03 | 1.053928e-03 | 6.353789e-14 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 4558 | 6.411814e-25 | 8.926255e-26 | 8.253532e-20 | 4.507500e-23 | 7.205970e-31 | 1.726429e-27 | 1.250789e-20 | 2.686227e-33 | 1.311215e-24 | 5.897132e-36 | 1.697311e-52 | 1.384834e-128 | 2.900642e-23 | 1.000000e+00 | 5.313257e-22 | 8.765537e-23 | 5.307327e-32 |
| 4559 | 2.086323e-06 | 3.369549e-17 | 3.327671e-04 | 5.864553e-07 | 1.086238e-07 | 1.224153e-08 | 5.838407e-07 | 5.278301e-09 | 2.777906e-07 | 2.360489e-10 | 2.657855e-17 | 2.734170e-98 | 1.430603e-05 | 9.996468e-01 | 2.142207e-06 | 3.534102e-07 | 3.758060e-10 |
| 4560 | 6.734257e-06 | 5.110178e-09 | 8.428272e-01 | 6.674259e-04 | 1.263140e-10 | 2.121082e-15 | 1.312298e-04 | 4.909820e-13 | 1.038616e-05 | 1.172244e-14 | 1.265739e-29 | 3.324217e-159 | 4.182381e-04 | 1.496179e-01 | 5.425747e-03 | 8.951116e-04 | 2.495943e-11 |
| 4561 | 7.538313e-09 | 4.016624e-17 | 4.030476e-04 | 4.667848e-02 | 3.162194e-10 | 8.476129e-09 | 3.475742e-06 | 2.490872e-12 | 2.172408e-08 | 8.886891e-13 | 8.700514e-18 | 2.442650e-111 | 9.528133e-01 | 9.868672e-05 | 2.594641e-06 | 4.280504e-07 | 4.469339e-10 |
| 4562 | 1.310695e-07 | 9.317343e-11 | 6.886927e-05 | 2.338925e-07 | 2.904712e-08 | 8.041276e-09 | 2.240113e-10 | 1.578051e-08 | 8.162015e-09 | 8.659953e-10 | 3.571252e-16 | 7.665526e-60 | 3.887733e-07 | 9.999298e-01 | 4.433498e-07 | 7.314155e-08 | 9.792083e-09 |
4563 rows × 17 columns
確率最大のcell typeを各細胞に割り当てる。
bdata.obs['CellAssign_prediction'] = predictions.idxmax(axis=1).values
UMAPで見ると以下のようにcell typeが分布している。
sc.pl.umap(
bdata,
color=['leiden_scVI', 'CellAssign_prediction'],
frameon=False,
ncols=2
)
マーカー遺伝子を使ったアノテーションがあんまりうまくいってないっぽいので、今度は「リファレンスデータからのラベル転移」でやってみる。
高精度なアノテーションが付与されたシングルセル解析データが大量に蓄積されているアトラスとして、マウスの場合はたとえばTabula Muris, Tabula Muris Senisがある。
ただどちらにも欲しい組織(網膜発生プロセス)が欲しい解像度で収載されていない。
ここでは、別の論文のデータを使う。網膜神経系発生について、胚から産後まで10個のタイムポイントで大量の細胞(10万以上)を対象にシングルセル解析したデータ。
Clark, Brian S., et al. "Single-cell RNA-seq analysis of retinal development identifies NFI factors as regulating mitotic exit and late-born cell specification." Neuron 102.6 (2019): 1111-1126. https://doi.org/10.1016/j.neuron.2019.04.010
この論文の素晴らしいところは、すべての細胞の遺伝子発現カウントデータとcell typeのアノテーションを公開してくれているところ。
マウス網膜神経発生アトラスといって差し支えないこのデータを活用して、ここまで扱ってきたデータとの統合解析を実行し、アトラスに付与されたcell typeアノテーションをラベル転移(Label transfer)する。
ラベル転移の半教師あり学習には scANVI を使う。確率モデルの詳細は以下の論文を参照。
Xu, Chenling, et al. "Probabilistic harmonization and annotation of single‐cell transcriptomics data with deep generative models." Molecular systems biology 17.1 (2021): e9620. https://doi.org/10.15252/msb.20209620
scANVIはscvi-toolsにモデルのひとつとして実装されている。
めっちゃデカいファイルで、GPU使わない場合は学習にも時間かかるので、データ読み込みと学習については講習では全部スキップ。学習されたモデルからラベルを予測した結果だけcsvファイルで配布する。
リファレンスデータの読み込み。マトリックス、遺伝子メタデータ、細胞メタデータがそれぞれ個別にファイルとして置いてあるので、それぞれ読み込んでひとつのanndataに統合する。
%%script true
mtx = sc.read_mtx('./ref_data/10x_mouse_retina_development.mtx')
refdata = mtx.transpose()
refdata.obs = pd.read_csv('./ref_data/10x_Mouse_retina_pData_umap2_CellType_annot_w_horiz.csv')
refdata.obs = refdata.obs.set_index('barcode')
refdata.var = pd.read_csv('./ref_data/10x_mouse_retina_development_feature.csv')
refdata.var = refdata.var.set_index('gene_short_name')
リファレンスデータのUMAP座標は計算済みのやつが細胞メタデータに記載されているので、scanpy.pl.umapで自動で読み込めるように、obsmの中に入れておく。このリファレンスデータの全体像がさくっと確認できる。
%%script true
refdata.obsm['X_umap'] = refdata.obs[['umap_coord1', 'umap_coord2']].values
sc.pl.umap(refdata,
color='umap2_CellType',
frameon=False)
クエリデータ(ここまで扱ってきたデータ)を カウントデータ として用意する。scANVIはscVIと同様、カウントデータを学習する確率モデルであるため。別レイヤーに取っておいたカウントのデータから新しいanndataを作る。
リファレンスデータは 'sample' のラベルのところに別々のタイムポイントからとったサンプルのラベルが記載されている。これをバッチと捉えてバッチ補正を行いたいので、クエリ側のバッチラベル("E2", "F2")もここに記述しておく。
%%script true
from anndata import AnnData
query = AnnData(X=adata.layers['counts'], obs=adata.obs, var=adata.var)
query.obs['sample'] = query.obs['batch']
リファレンスとクエリのデータをひとつのanndataにまとめる。
%%script true
refdata.var_names_make_unique()
alldata = refdata.concatenate(query)
前半でやったのと同じような、シングルセル解析の前処理をデータ全体に対して実行する。
%%script true
alldata.layers['counts'] = alldata.X.copy()
sc.pp.normalize_per_cell(alldata, counts_per_cell_after=1e4)
sc.pp.log1p(alldata)
alldata.raw = alldata
%%script true
sc.pp.highly_variable_genes(alldata, flavor='seurat_v3', n_top_genes=2000, layer='counts', batch_key='sample')
scVIを利用してデータ全体のバッチ補正を行う統合モデルを学習する。CPUで学習する場合は時間かかるので昼食でもとりにいく。帰ってきたら学習終わってる、くらいの計算時間で済むはず。
%%script true
scvi.model.SCVI.setup_anndata(alldata, layer='counts', batch_key='sample')
%%script true
integration_model = scvi.model.SCVI(alldata)
integration_model.train()
%%script true
integration_model.save('./models/integration_model')
%%script true
integration_model = scvi.model.SCVI.load('./models/integration_model', alldata)
リファレンス側のデータでcell typeが記述されているカラム('umap2_CellType')が、クエリ側のデータでは現在 NaN になっているので、ここを全部 'Unknown' にセットしておく。
%%script true
alldata.obs['umap2_CellType'] = alldata.obs['umap2_CellType'].cat.add_categories('Unknown')
alldata.obs = alldata.obs.fillna({'umap2_CellType':'Unknown'})
統合モデル、統合データ、未知ラベルのカラムと名前を指定してscANVIモデルをセットアップする。
%%script true
label_model = scvi.model.SCANVI.from_scvi_model(integration_model,
adata=alldata,
unlabeled_category='Unknown',
labels_key='umap2_CellType')
ラベル転移の学習を実行する。ここはそれほど計算時間かからない(ようにパラメータを設定しているが精度はじゅうぶん出る)
%%script true
label_model.train(n_samples_per_label=100)
label_model.save('./models/label_transfer')
%%script true
label_model = scvi.model.SCANVI.load('./models/label_transfer/', alldata)
学習されたモデルでcell typeの予測を実行する。
%%script true
alldata.obs['predicted_celltype'] = label_model.predict(alldata)
この予測結果全体からクエリ側に対応するデータだけ抜き出して、以下のようにcsvとして保存したのが配布しているデータ。
%%script true
predictions = alldata.obs[(alldata.obs['sample'] == 'E2') | (alldata.obs['sample'] == 'F2')]['predicted_celltype']
predictions.index = predictions.index.str[:-2]
predictions.to_csv('./models/label_transfer/celltype_predictions.csv')
csvファイルをロード。
predictions = pd.read_csv('./models/label_transfer/celltype_predictions.csv', index_col=0)
adataに予測されたラベルを格納。
adata.obs['predicted_celltype'] = predictions.loc[adata.obs.index, 'predicted_celltype']
予測されたラベルの細胞数カウントは以下。それっぽいラベルが並んでいる。
adata.obs['predicted_celltype'].value_counts()
Retinal Ganglion Cells 1198 Late RPCs 1118 Neurogenic Cells 635 Amacrine Cells 517 Early RPCs 489 Horizontal Cells 263 Photoreceptor Precursors 218 Cones 102 Bipolar Cells 16 Rods 7 Name: predicted_celltype, dtype: int64
以下が予測結果で色分けしたUMAP。
RPCs は Retinal Progenitor Cells の略。
RPCs からNeurogenicにつながって、Photoreceptor/Conesに分岐していく流れと、Amacrine/Horizontalに分岐していく流れと、Retinal ganglion cellsに分岐する3つの流れに分かれていることがわかった。
sc.pl.umap(adata,
color=['leiden_scVI', 'predicted_celltype'],
frameon=False, alpha=0.5,
ncols=2)
Leidenクラスタに名前をつけたいので、それぞれのクラスタに所属している細胞がどのcell typeラベルを持っているのか、数をクロス集計してみる。
pd.crosstab(adata.obs['leiden_scVI'], adata.obs['predicted_celltype'])
| predicted_celltype | Amacrine Cells | Bipolar Cells | Cones | Early RPCs | Horizontal Cells | Late RPCs | Neurogenic Cells | Photoreceptor Precursors | Retinal Ganglion Cells | Rods |
|---|---|---|---|---|---|---|---|---|---|---|
| leiden_scVI | ||||||||||
| 0 | 0 | 10 | 2 | 177 | 3 | 536 | 63 | 0 | 0 | 0 |
| 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 600 | 0 |
| 2 | 359 | 0 | 0 | 0 | 122 | 0 | 9 | 16 | 0 | 0 |
| 3 | 1 | 1 | 1 | 100 | 0 | 382 | 7 | 0 | 0 | 0 |
| 4 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 488 | 0 |
| 5 | 0 | 3 | 0 | 205 | 0 | 198 | 15 | 0 | 0 | 0 |
| 6 | 0 | 0 | 0 | 2 | 1 | 2 | 320 | 0 | 0 | 0 |
| 7 | 0 | 0 | 99 | 2 | 0 | 0 | 0 | 151 | 0 | 7 |
| 8 | 0 | 1 | 0 | 2 | 86 | 0 | 36 | 10 | 109 | 0 |
| 9 | 0 | 0 | 0 | 1 | 8 | 0 | 184 | 41 | 1 | 0 |
| 10 | 156 | 1 | 0 | 0 | 42 | 0 | 0 | 0 | 0 | 0 |
このカウントを元に、多数決で以下のような適当な名前をつけた。
leiden_to_celltypes = {
'0':'Late-RPCs-1', '1':'RGCs-3',
'2':'A/H-1', '3':'Late-RPCs-2', '4':'RGCs-2',
'5':'Early-RPCs', '6':'Neurogenic-1',
'7':'Photoreceptors',
'8':'RGCs-1',
'9':'Neurogenic-2',
'10':'A/H-2'}
adata.obs['leiden_celltypes'] = [leiden_to_celltypes[v] for v in adata.obs['leiden_scVI'].values]
名前がついたLeidenクラスタを可視化してみる。
with plt.rc_context({"figure.figsize": (9, 9)}):
sc.pl.umap(adata,
color='leiden_celltypes',
legend_fontsize=10,
frameon=False, alpha=0.5,
legend_loc='on data')
scVeloを使うためには、spliced/unsplicedの二種類のカウントデータが必要となる。
なので事前に velocyto を使って、Cell Rangerが生成したゲノムマッピングのBAMファイルから、spliced/unsplicedのカウントデータを作っておく。
バッチごとに実行して、それぞれのバッチに属する細胞のspliced/unsplicedカウントを格納した loom データを作る。
計算に時間がかかるので、講習ではすでにvelocyto計算済みの結果を配布している。
velocytoは以下のコマンドで実行する。
velocyto run10x \
-m /path/to/Cell_Ranger_References/mm10_rmsk.gtf \
/path/to/RetinalBatchE2 \
/path/to/Cell_Ranger_References/refdata-gex-mm10-2020-A/genes/genes.gtf
velocyto run10x \
-m /path/to/Cell_Ranger_References/mm10_rmsk.gtf \
/path/to/RetinalBatchF2 \
/path/to/Cell_Ranger_References/refdata-gex-mm10-2020-A/genes/genes.gtf
velocytoでバッチごとに推定したloomファイルを統合して書き出す。ここの統合処理も時間かかるので、統合後のデータを配布。
%%script true
import loompy
loompy.combine(['./data/RetinalBatchE2/velocyto/RetinalBatchE2.loom',
'./data/RetinalBatchF2/velocyto/RetinalBatchF2.loom'],
output_file='./data/retinal_velo.loom')
import scvelo as scv
scv.set_figure_params()
spliced/unsplicedのカウントデータを遺伝子発現テーブルと統合する。細胞のバーコードで紐づけるので、長さの指定(インデックスのどの部分がバーコードか)が必要になることに注意。
ldata = scv.read('./data/retinal_velo.loom', cache=True)
# id_lengthは細胞を区別するバーコード配列の長さ。事前に確認しておく
# また、scVIのmerge関数は勝手にこれまでのクラスタの色指定をリセットしちゃうので一時退避させて再設定しておく
# 'uns'(unstructured)は特に構造の決まってない雑多なメタデータを格納しておく場所。
cluster_colors = adata.uns['leiden_celltypes_colors']
adata = scv.utils.merge(adata, ldata, id_length=16)
adata.uns['leiden_celltypes_colors'] = cluster_colors
... reading from cache file cache/data-retinal_velo.h5ad
scVeloの関数で、遺伝子ごとのspliced/unsplicedカウントの比率を表示できる。実験プラットフォームにもよるが、だいたいunsplicedが25%程度らしい。クラスタごとにも表示。極端にunsplicedがとれてないクラスタがあるかどうかチェックする。
scv.pl.proportions(adata, groupby='leiden_celltypes')
速度計算を実行する。近傍グラフの構成から。
# scVI補完の潜在空間上で近傍グラフ構成、一次・二次モーメント計算
scv.pp.moments(adata, use_rep='X_scVI')
WARNING: Did not normalize X as it looks processed already. To enforce normalization, set `enforce=True`.
Normalized count data: spliced, unspliced.
computing moments based on connectivities
finished (0:00:03) --> added
'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
Splicing kineticsモデルのパラメータ推論。
scv.tl.velocity(adata, mode='stochastic')
computing velocities
finished (0:00:04) --> added
'velocity', velocity vectors for each individual cell (adata.layers)
速度ベクトルを元に、細胞から細胞への遷移確率を計算。
scv.tl.velocity_graph(adata)
computing velocity graph (using 1/10 cores)
0%| | 0/4563 [00:00<?, ?cells/s]
finished (0:00:15) --> added
'velocity_graph', sparse matrix with cosine correlations (adata.uns)
速度ベクトルから構成した「流れ」をプロットしてみる。
scv.pl.velocity_embedding_stream(adata,
basis='X_umap', color='leiden_celltypes',
legend_fontsize=9,
smooth=0.8, min_mass=4)
computing velocity embedding
finished (0:00:00) --> added
'velocity_umap', embedded velocity vectors (adata.obsm)
Progenitorの集団からはじまって、わかりやすい三分岐。
RNA速度を計算することによって、トランスクリプトームパターンの幾何学的な近接性だけで推定された 'Pseudo time' ではなく、速度情報をちゃんと織り込んだ細胞の 'Latent time' を計算できる。
scv.tl.recover_dynamics(adata, n_jobs=1)
scv.tl.recover_latent_time(
adata,
root_key='initial_states_probs',
end_key='terminal_states_probs')
recovering dynamics (using 1/10 cores)
0%| | 0/798 [00:00<?, ?gene/s]
finished (0:06:01) --> added
'fit_pars', fitted parameters for splicing dynamics (adata.var)
computing terminal states
identified 0 region of root cells and 1 region of end points .
finished (0:00:00) --> added
'root_cells', root cells of Markov diffusion process (adata.obs)
'end_points', end points of Markov diffusion process (adata.obs)
WARNING: No root cells detected. Consider specifying root cells to improve latent time prediction.
computing latent time using root_cells as prior
finished (0:00:02) --> added
'latent_time', shared time (adata.obs)
scv.pl.scatter(
adata,
color=['leiden_celltypes', 'latent_time'],
fontsize=16,
cmap='viridis',
perc=[2, 98],
colorbar=True,
rescale_color=[0, 1],
title=['clusters', 'latent time'])
Progenitorの部分がごちゃごちゃしてわかりにくいので、取り出して個別に解析してみる。
どうもグルグルとまわってるみたいなので、細胞周期のマーカー遺伝子発現をチェック。
prog = adata[adata.obs['leiden_celltypes'].isin(['Early-RPCs', 'Late-RPCs-1', 'Late-RPCs-2'])].copy()
sc.pp.neighbors(prog, use_rep="X_scVI", n_neighbors=50)
sc.tl.umap(prog, min_dist=1.0)
sc.pl.umap(prog,
ncols=2,
color=['Mcm6', 'Esco2', 'Top2a', 'Aurka', 'Cenpa'],
cmap='viridis',
frameon=False)
computing neighbors
finished: added to `.uns['neighbors']`
`.obsp['distances']`, distances for each pair of neighbors
`.obsp['connectivities']`, weighted adjacency matrix (0:00:03)
computing UMAP
finished: added
'X_umap', UMAP coordinates (adata.obsm) (0:00:03)
というわけで、このグルグルは細胞周期を反映しているのだった。
この中だけで速度を描くとよりわかりやすい。
scv.pp.moments(prog, use_rep='X_scVI')
scv.tl.velocity(prog, mode='stochastic')
scv.tl.velocity_graph(prog)
scv.pl.velocity_embedding_stream(prog,
basis='X_umap', color='leiden_celltypes',
smooth=0.8, min_mass=4)
computing moments based on connectivities
finished (0:00:01) --> added
'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
computing velocities
finished (0:00:01) --> added
'velocity', velocity vectors for each individual cell (adata.layers)
computing velocity graph (using 1/10 cores)
0%| | 0/1704 [00:00<?, ?cells/s]
finished (0:00:02) --> added
'velocity_graph', sparse matrix with cosine correlations (adata.uns)
computing velocity embedding
finished (0:00:00) --> added
'velocity_umap', embedded velocity vectors (adata.obsm)
「膨らみ方」(splicing kinetics)は遺伝子ごとに異なる。scVeloのscatter関数では遺伝子ごとに、具体的にどのようなダイナミクスが推定されたのかプロットしてくれる。それぞれの遺伝子がそれぞれ異なるタイミング、異なるパターンで、induction => steady-state => repressionのパスを巡っている。
scv.pl.scatter(prog,
['Mcm6', 'Esco2', 'Top2a', 'Aurka', 'Cenpa'],
ncols=2,
color='leiden_celltypes')
細かく見ると印象が違ってくることもあるので、流れの全体像として粗視化されたパターンだけでなく、それぞれの速度をつぶさに見ていくことも大事。
scv.pl.velocity_embedding(adata, figsize=(9, 9),
basis='X_umap', color='leiden_celltypes',
scale=0.5)
Trajectoryを推定するツールは数多あるが、幾何学的な近接性だけではなくRNA速度情報も考慮に入れて道筋を推定するツールはあまり多くない。
ここでは、幾何学的近接性をベースにしたTrajectory解析のロバスト性と、RNA速度による方向情報を組み合わせることによって、方向づけされた、確率的な状態変化の軌跡を推定するツール、 CellRank を使ってみる。
Lange, Marius, et al. "CellRank for directed single-cell fate mapping." Nature methods 19.2 (2022): 159-170. https://doi.org/10.1038/s41592-021-01346-6
近傍グラフ上で細胞それぞれを状態としたマルコフ連鎖を考え、速度ベクトルと遺伝子発現パターンの類似性のふたつを加味して、細胞間の遷移確率を推定、それにより細胞系譜の初期・中間・終端の細胞集団を特定して、細胞運命決定の確率的な性質とRNA速度推定の不確かさを織り込んだ、細胞の「運命確率」を計算できる。また、細胞系譜を駆動する遺伝子を特定し、その系譜に特異的な遺伝子の発現トレンドを明らかにすることができる。
CellRankの計算には、事前にscVeloなどのツールを利用して速度ベクトルを計算しておくことが必要。
import cellrank as cr
import numpy as np
cr.settings.verbosity = 2
細胞系譜の「終端」の状態に該当する細胞集団を特定する。
cr.tl.terminal_states(adata,
cluster_key='leiden_celltypes',
weight_connectivities=0.5,
n_states=4,
mode='deterministic',
force_recompute=True)
Computing transition matrix based on logits using `'deterministic'` mode Estimating `softmax_scale` using `'deterministic'` mode
0%| | 0/4563 [00:00<?, ?cell/s]
Setting `softmax_scale=11.4320`
0%| | 0/4563 [00:00<?, ?cell/s]
Finish (0:00:19)
Using a connectivity kernel with weight `0.5`
Computing transition matrix based on `adata.obsp['connectivities']`
Finish (0:00:00)
Computing Schur decomposition
Adding `adata.uns['eigendecomposition_fwd']`
`.schur_vectors`
`.schur_matrix`
`.eigendecomposition`
Finish (0:00:00)
Computing `4` macrostates
Adding `.macrostates`
`.macrostates_memberships`
`.coarse_T`
`.coarse_initial_distribution
`.coarse_stationary_distribution`
`.schur_vectors`
`.schur_matrix`
`.eigendecomposition`
Finish (0:00:00)
Adding `adata.obs['terminal_states']`
`adata.obs['terminal_states_probs']`
`.terminal_states`
`.terminal_states_probabilities`
`.terminal_states_memberships
Finish`
cr.pl.terminal_states(adata)
細胞系譜の「初期」状態に該当する細胞集団を特定する。
cr.tl.initial_states(adata,
cluster_key='leiden_celltypes',
weight_connectivities=0.5,
force_recompute=True)
cr.pl.initial_states(adata, discrete=True)
Computing transition matrix based on logits using `'deterministic'` mode Estimating `softmax_scale` using `'deterministic'` mode
0%| | 0/4563 [00:00<?, ?cell/s]
Setting `softmax_scale=11.4320`
0%| | 0/4563 [00:00<?, ?cell/s]
Finish (0:00:15)
Using a connectivity kernel with weight `0.5`
Computing transition matrix based on `adata.obsp['connectivities']`
Finish (0:00:00)
Computing eigendecomposition of the transition matrix
Adding `adata.uns['eigendecomposition_bwd']`
`.eigendecomposition`
Finish (0:00:00)
Computing Schur decomposition
Adding `adata.uns['eigendecomposition_bwd']`
`.schur_vectors`
`.schur_matrix`
`.eigendecomposition`
Finish (0:00:00)
Computing `2` macrostates
Adding `.macrostates`
`.macrostates_memberships`
`.coarse_T`
`.coarse_initial_distribution
`.coarse_stationary_distribution`
`.schur_vectors`
`.schur_matrix`
`.eigendecomposition`
Finish (0:00:00)
Adding `adata.obs['initial_states']`
`adata.obs['initial_states_probs']`
`.terminal_states`
`.terminal_states_probabilities`
`.terminal_states_memberships
Finish`
推定された終端状態に向けて、それぞれの細胞が「どの終端にどれだけの確率で向かうか」を推定する。
cr.tl.lineages(adata)
cr.pl.lineages(adata,
ncols=2,
same_plot=False)
Computing absorption probabilities
0%| | 0/4 [00:00<?, ?/s]
Adding `adata.obsm['to_terminal_states']`
`.absorption_probabilities`
Finish (0:00:00)
細胞系譜を駆動する遺伝子を見つける。
これは単純にピアソン相関でランキングしているらしい。それぞれの終端に向けた運命確率が計算済みなので、細胞ごとに運命確率と遺伝子発現量の相関を計算して、それが高い遺伝子が細胞系譜を駆動している、という仮定。
結果は全遺伝子について、それぞれの終端ごとの相関係数とピアソン相関のp-valueなどが入ったデータフレームとして出力される。
cr.tl.lineage_drivers(adata,
cluster_key='leiden_celltypes')
Adding `adata.varm['terminal_lineage_drivers']`
`.lineage_drivers`
Finish (0:00:00)
| Photoreceptors_corr | Photoreceptors_pval | Photoreceptors_qval | Photoreceptors_ci_low | Photoreceptors_ci_high | A/H-1_corr | A/H-1_pval | A/H-1_qval | A/H-1_ci_low | A/H-1_ci_high | RGCs-2_corr | RGCs-2_pval | RGCs-2_qval | RGCs-2_ci_low | RGCs-2_ci_high | Late-RPCs-2_corr | Late-RPCs-2_pval | Late-RPCs-2_qval | Late-RPCs-2_ci_low | Late-RPCs-2_ci_high | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Gngt2 | 0.745112 | 0.000000e+00 | 0.000000e+00 | 0.731920 | 0.757746 | -0.375494 | 1.361511e-156 | 4.033704e-154 | -0.400151 | -0.350295 | -0.111782 | 3.453654e-14 | 7.790882e-13 | -0.140343 | -0.083035 | -0.297182 | 4.068856e-95 | 1.208151e-93 | -0.323409 | -0.270498 |
| Crx | 0.712121 | 0.000000e+00 | 0.000000e+00 | 0.697517 | 0.726133 | -0.349282 | 7.283963e-134 | 1.618496e-131 | -0.374502 | -0.323544 | -0.109482 | 1.146391e-13 | 2.473089e-12 | -0.138060 | -0.080722 | -0.292433 | 5.574067e-92 | 1.594709e-90 | -0.318745 | -0.265671 |
| Thrb | 0.696499 | 0.000000e+00 | 0.000000e+00 | 0.681250 | 0.711143 | -0.353160 | 4.395616e-137 | 1.010385e-134 | -0.378300 | -0.327500 | -0.077313 | 1.684023e-07 | 1.925506e-06 | -0.106091 | -0.048405 | -0.292105 | 9.133293e-92 | 2.601818e-90 | -0.318422 | -0.265337 |
| Neurod4 | 0.677728 | 0.000000e+00 | 0.000000e+00 | 0.661725 | 0.693115 | -0.243792 | 2.397837e-63 | 1.402104e-61 | -0.270892 | -0.216306 | -0.128552 | 2.566051e-18 | 7.739952e-17 | -0.156983 | -0.099908 | -0.356125 | 1.409289e-139 | 7.171236e-138 | -0.381202 | -0.330524 |
| Gnb3 | 0.631526 | 0.000000e+00 | 0.000000e+00 | 0.613757 | 0.648656 | -0.306965 | 8.847064e-102 | 1.123324e-99 | -0.333015 | -0.280446 | -0.095507 | 9.861265e-11 | 1.625098e-09 | -0.124179 | -0.066675 | -0.263217 | 4.891180e-74 | 1.088635e-72 | -0.290019 | -0.236004 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| Sox11 | -0.331815 | 5.856069e-120 | 1.419511e-117 | -0.357391 | -0.305742 | 0.589021 | 0.000000e+00 | 0.000000e+00 | 0.569742 | 0.607652 | 0.219928 | 1.672865e-51 | 2.186534e-49 | 0.192137 | 0.247365 | -0.412515 | 8.201189e-193 | 6.469719e-191 | -0.436309 | -0.388144 |
| Ptma | -0.336005 | 3.302656e-123 | 8.467501e-121 | -0.361497 | -0.310011 | 0.299929 | 5.834488e-97 | 6.648325e-95 | 0.273291 | 0.326107 | 0.213032 | 2.450786e-48 | 2.722823e-46 | 0.185160 | 0.240561 | -0.101739 | 5.421013e-12 | 2.429343e-11 | -0.130371 | -0.072938 |
| Eno1 | -0.356491 | 6.897951e-140 | 2.090079e-137 | -0.381561 | -0.330898 | -0.147829 | 8.619659e-24 | 1.427544e-22 | -0.176090 | -0.119325 | -0.017537 | 2.362626e-01 | 4.149423e-01 | -0.046530 | 0.011485 | 0.527944 | 0.000000e+00 | 0.000000e+00 | 0.506689 | 0.548557 |
| Dkk3 | -0.399434 | 1.961110e-179 | 8.434040e-177 | -0.423541 | -0.374761 | 0.052400 | 3.975221e-04 | 1.357174e-03 | 0.023420 | 0.081293 | -0.247073 | 4.378924e-65 | 9.121846e-63 | -0.274125 | -0.219632 | 0.502244 | 1.921847e-304 | 3.713342e-302 | 0.480226 | 0.523629 |
| Pax6 | -0.452822 | 2.181607e-238 | 1.454259e-235 | -0.475590 | -0.429448 | 0.502436 | 1.007049e-304 | 1.678248e-301 | 0.480423 | 0.523815 | 0.032143 | 2.990866e-02 | 9.101626e-02 | 0.003130 | 0.061102 | -0.084340 | 1.137229e-08 | 4.078972e-08 | -0.113080 | -0.055460 |
13332 rows × 20 columns
終端ごとにそれぞれの発現量をプロットしてみる。
cr.pl.lineage_drivers(adata,
lineage='Photoreceptors', n_genes=4, ncols=4)
cr.pl.lineage_drivers(adata,
lineage='A/H-1', n_genes=4, ncols=4)
cr.pl.lineage_drivers(adata,
lineage='RGCs-2', n_genes=4, ncols=4)
細胞が辿る終端に向けた道筋の上で、遺伝子発現がどのようなトレンドで推移しているのかを計算し、プロットできる。
単純にscVeloで計算したlatent timeに対する全細胞の平均的な発現量を見るのではなく、それぞれの運命への所属確率で発現量を重み付けして表現している。
model = cr.ul.models.GAM(adata)
cr.pl.gene_trends(
adata,
model=model,
data_key='X',
genes=['Gngt2', 'Tfap2b', 'Pou4f1'],
ncols=3,
time_key="latent_time",
same_plot=True,
hide_cells=True,
figsize=(15, 4),
n_test_points=200)
Computing trends using `1` core(s)
0%| | 0/3 [00:00<?, ?gene/s]
Finish (0:00:00) Plotting trends
また、lineage driverとして特定された複数の遺伝子の発現トレンドの推移をヒートマップとしてプロットすることもできる。
cr.pl.heatmap(
adata,
model,
genes=adata.varm['terminal_lineage_drivers']['A/H-1_corr'].sort_values(ascending=False).index[:100],
lineages='A/H-1')
Computing trends using `1` core(s)
0%| | 0/100 [00:00<?, ?gene/s]
Finish (0:00:05)